#include <iostream>
#include <cmath>


const double PI = (std::atan(1.0) * 4.0);

double source(double x, double y);

double u_b(double x, double y);

double source(double x, double y)
{
    return 2.0 * PI * PI * std::sin(PI * x) * std::sin(PI * y);
};

double u_b(double x, double y)
{
    return std::sin(PI * x) * std::sin(PI * y);
};

int main(int argc, char *argv[])
{
    double x0 = 0.0;
    double y0 = 0.0;
    double x1 = 1.0;
    double y1 = 1.0;

    int N = 20;
    double h = (x1 - x0) / (N - 1);
    int m = N - 2;
    int n = m * m;

    double A[n][n];
    double f[n];

    for (int i = 0; i < n; i++)
    {
	for (int j = 0; j < n; j++)
	    A[i][j] = 0;
	f[i] = 0;
    }

    for (int i = 0; i < m; i++)
	for (int j = 0; j < m; j++)
	{
	    int idx = i * m + j;
	    double lx = (j + 1) * h;
	    double ly = (i + 1) * h;
	    A[idx][idx] = 4;
	    f[idx] = source(lx, ly);
	    if (j != 0)
		A[idx][idx - 1] = -1;
	    else
		f[idx] += u_b(lx - h, ly);
	    if (j != m - 1)
		A[idx][idx + 1] = -1;
	    else
		f[idx] += u_b(lx + h, ly);
	    if (i != 0)
		A[idx][idx - m] = -1;
	    else
		f[idx] += u_b(lx, ly - h);
	    if (i != m - 1)
		A[idx][idx + m] = -1;
	    else
		f[idx] += u_b(lx, ly + h);
	}
    std::cout << "A = [" << std::endl;
    for (int i = 0; i < n; i++)
    {
	for (int j = 0; j < n; j++)
	    std::cout << A[i][j] << "\t";
	std::cout << std::endl;
    }
    std::cout << "];" << std::endl;
    std::cout << "f = [" << std::endl;
    for (int i = 0; i < n; i++)
	std::cout << f[i] << std::endl;
    std::cout << "];" << std::endl;
    return 0;
};
